package jmars.mobile.rbscluster;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileWriter;
import java.io.IOException;
import java.io.InputStreamReader;
import java.util.ArrayList;
import java.util.Date;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

import org.apache.commons.math3.stat.StatUtils;

import weka.core.matrix.Matrix;
import weka.core.matrix.EigenvalueDecomposition;

/**
 * Construct pairwise affinity matrix A i,j = exp(-|xi-xj|/2sigma^2)
Construct degree matrix D = diag(d1,,,dn)
Compute Laplacian L=D-A(unormalised)
	Lsym=D^(-1/2)LD^(-1/2)(normalized)
Compute the first k eigen-vectors u1,,,uk of L
Let U in R_nxk contain the vectors u1,,,uk as columns
Let yi in Rk be the vector corresponding to the i-th row of U
Cluster the points (Vi) i=1,,,N into k clusters h1,,,,hk with K-means
 * @author Jan
 *
 */
public class SpectralKMeansCluster {

	//private static String path="C:\\Users\\ma01m\\workspace\\datainvestigation\\data\\";
	private static String path="E:\\BaiduYunDownload\\afterprocess\\";
	private static int samplingSkip = 1;
	private static String module="yuzhong";
	private static String loadFile=path+module+"\\rbs_cluster1080";//1054 1080 1130
	
	private static int clusterNum=100;
	private static Map<String,String[]> instances=new HashMap<String,String[]>();
	private static Map<String,double[]> trafficRank = new HashMap<String,double[]>();
	private static Map<String,double[]> normalizedTraffic = new HashMap<String,double[]>();
	private static List<String> instanceKeyList = new ArrayList<String>();
	private static Matrix M;
	private static Map<Integer,String[]> clusters=new HashMap<Integer,String[]>();
	private static Map<String,Integer> results= new HashMap<String,Integer>();

	static double[][] L,D;
	static FileWriter resultWriter;
	static Matrix Md,Ml,MdSqrtInverse,Mlsym,Mlrw;
	
	private static double affinityThreshold=0;
	private static String lFlag_RW="rw";
	private static String lFlag_SYM="sym";
	private static String lFlag=lFlag_SYM;
	
	private static String getFilename(){
		return path+module+"\\mat_cluster_data_1_of_"+samplingSkip+"_"+(affinityThreshold==0?"fullgraph":"knn")+"_"+sigmaS;
	}
	
	private static double sigma;
	private static String sigmaS;
	
	public static void main(String[] args) throws IOException {
		for(int t=10;t<=20;t+=5){
			sigma=(double)t/100;
			sigmaS="0."+t;
			System.out.println(new Date()+": starting with sigma of "+sigma+"*****************************");
			run();
			System.gc();
		}
		System.out.println(new Date()+": finished*****************************");
	}
	public static void run() throws IOException{

		// load instances
		System.out.println(new Date()+": start");
		loadInstances();
		System.out.println(new Date()+": instances loaded");
		
		/*
		System.out.println(affinity(0.05,sigma));
		System.out.println(affinity(0.15,sigma));
		System.out.println(affinity(0.25,sigma));
		System.out.println(affinity(0.35,sigma));
		System.out.println(affinity(0.45,sigma));
		System.out.println(affinity(0.55,sigma));
		System.out.println(affinity(0.65,sigma));
		System.out.println(affinity(0.75,sigma));
		System.out.println(affinity(0.85,sigma));
		System.out.println(affinity(0.95,sigma));
		System.out.println(affinity(1.05,sigma));
		System.out.println(affinity(1.15,sigma));
		System.out.println(affinity(1.25,sigma));
		System.out.println(affinity(1.35,sigma));
		System.out.println(affinity(1.45,sigma));
		System.out.println(affinity(1.56,sigma));
		System.exit(0);
		*/
//		testDistance();
//		System.exit(0);
		
		// build Laplacian matrix
		L=new double[instanceKeyList.size()][instanceKeyList.size()];
		//unnormalized Laplacian matrix
		for(int i=0;i<instanceKeyList.size();i++){
			for(int j=0;j<i;j++){
				L[i][j]=-affinity(distance(instances.get(instanceKeyList.get(i)),instances.get(instanceKeyList.get(j))),sigma);
				L[j][i]=L[i][j];
				L[i][i]=L[i][i]-L[i][j];
				L[j][j]=L[j][j]-L[j][i];
			}
		}
		D=new double[instanceKeyList.size()][instanceKeyList.size()];
		//Degree matrix
		for(int i=0;i<instanceKeyList.size();i++){
			D[i][i]=L[i][i];
		}
		System.out.println(new Date()+": unnormalized Laplacian matrix built");
		
		//build normalized Laplacian matrix
		Ml=new Matrix(L);
		Md=new Matrix(D);
		
		if(lFlag.equals(lFlag_SYM)){
			MdSqrtInverse = Md.sqrt().inverse();
			Mlsym=MdSqrtInverse.times(Ml).times(MdSqrtInverse);
			M=Mlsym;
		}else{
			Mlrw=Md.inverse().times(Ml);
			M=Mlrw;
		}
		System.out.println(new Date()+": normalized Laplacian matrix built");
		
		// solve eigen value & vector
		EigenvalueDecomposition E = new EigenvalueDecomposition(M);
		System.out.println(new Date()+": eigen values & vetors decomposited");
		
		// dimension reduction & output for matlab
		resultWriter = new FileWriter(getFilename()+"_eigenvalues");
		int[] order = new int[E.getRealEigenvalues().length];
		int temp;
		for(int t=0;t<order.length;t++)
			order[t]=t;
		for(int i=0;i<order.length;i++){
			for(int j=0;j<order.length-1-i;j++){
				if(E.getRealEigenvalues()[order[j]]>E.getRealEigenvalues()[order[j+1]]){
					temp=order[j];order[j]=order[j+1];order[j+1]=temp;
				}
			}
		}
		resultWriter.write("%eigenValues:"+E.getRealEigenvalues().length+"\r\n");
		for(int t=0;t<order.length;t++){
			if(Math.abs(E.getRealEigenvalues()[order[t]])<1E-4)
				resultWriter.write("0\r\n");
			else
				resultWriter.write(E.getRealEigenvalues()[order[t]]+"\r\n");
		}
		resultWriter.close();
		System.out.println(new Date()+": eigen values output finished");
//		for(int t=0;t<clusterNum;t++){
//			resultWriter = new FileWriter(getFilename()+"_eigenvectors_"+t);
//			resultWriter.write("%eigenvalue:"+E.getRealEigenvalues()[order[t]]+"\r\n");
//			for(int i=0;i<instanceKeyList.size();i++){
//				resultWriter.write(ev.get(i, t)+"\r\n");
//			}
//			resultWriter.close();
//		}
//		System.out.println(new Date()+": eigen vectors output finished");
		resultWriter = new FileWriter(getFilename());
		resultWriter.write("%@RELATION rbs\r\n");
		for(int t=0;t<clusterNum;t++)
			resultWriter.write("%@Attribute ev"+t+" numeric\r\n");
		resultWriter.write("%@data\r\n");
		Matrix ev = E.getV();
		double tempDouble;
		for(int i=0;i<instanceKeyList.size();i++){
			if(lFlag.equals(lFlag_SYM)){
				//normalizing the rows to norm 1, when we use Lsym
				tempDouble=0;
				for(int j=0;j<clusterNum;j++){
					tempDouble+=Math.pow(ev.getArray()[i][order[j]],2);
				}
				tempDouble=Math.sqrt(tempDouble);
				for(int j=0;j<clusterNum;j++){
					resultWriter.write(ev.getArray()[i][order[j]]/tempDouble+((j!=clusterNum-1)?",":"\r\n"));
				}
			}else{
				for(int j=0;j<clusterNum;j++){
					resultWriter.write(ev.getArray()[i][order[j]]+((j!=clusterNum-1)?",":"\r\n"));
				}
			}
		}
		resultWriter.close();
		System.out.println(new Date()+": clustering data output finished");
		// clustering
//		
//		double[] x={1,2,3,5,8};
//		double[] y={1,2,3,4,8};
//		System.out.println(pearsonRho(x,y));
//		System.out.println(spearmanRho(x,y));
		
	}
	


	public static void cluster(){
		Double minDistance=null,d=null;
		for(String key:instances.keySet()){
			if(clusters.size()<clusterNum){
				results.put(key, clusters.size());
				clusters.put(clusters.size(), instances.get(key));
				continue;
			}
			d=0D;//calculate d
			minDistance=null;
			for(int t=0;t<clusters.size();t++){
				if(minDistance==null) {
					minDistance = d;
					results.put(key, t);
				}else if(minDistance>d){
					minDistance = d;
					results.put(key, t);
				}
			}
		}
		d=calClusterCenter();
	}
	
	public static double calClusterCenter(){
		Double d=0D;
		for(@SuppressWarnings("unused") String key:instances.keySet()){
			
		}
		return d;
	}
	
	
	
	
	
	private static void loadInstances(){
		try{
			String row;
			int temp;
			String[] cols;
			String lastTraffic;
			Integer lastIdx;
			File file= new File(loadFile);
			resultWriter = new FileWriter(loadFile+"_"+samplingSkip);
			FileInputStream ins = new FileInputStream(file);
			BufferedReader reader = new BufferedReader(new InputStreamReader(ins));
			instances.clear();
			instanceKeyList.clear();
			int count=0;
			while((row=reader.readLine())!=null){
				if(row.startsWith("%")||row.startsWith("@"))continue;
				if(row.trim().equals(""))continue;
				count++;
				if(count%samplingSkip!=0)continue;
				resultWriter.write(row+"\r\n");
				cols=row.split(",");
				instances.put(cols[0],cols);
				instanceKeyList.add(cols[0]);


				double[] traff=new double[24];
				
				int[] order = {2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25};
				for(int i=0;i<24;i++){
					for(int j=0;j<23-i;j++){
						if(Double.parseDouble(cols[order[j]])>Double.parseDouble(cols[order[j+1]])){
							temp=order[j];order[j]=order[j+1];order[j+1]=temp;
						}
					}
					traff[i]=Double.parseDouble(cols[i+2]);
				}
				
				double mean=StatUtils.mean(traff);
				double std=Math.sqrt(StatUtils.variance(traff));
				for(int t=0;t<24;t++){
					traff[t]=(traff[t]-mean)/std;
				}
				
				normalizedTraffic.put(cols[0], traff);
//				for(int t=0;t<24;t++) System.out.print(cols[t+2]+"  ");
//				System.out.println();
//				for(int t=0;t<24;t++) System.out.print(order[t]+"  ");
//				System.out.println();
//				for(int t=0;t<24;t++) System.out.print(cols[order[t]]+"  ");
//				System.out.println();
				
				double[] rank=new double[24];
				lastTraffic=null;
				lastIdx=null;
				for(int t=0;t<24;t++){
					if(lastTraffic==null){
						lastTraffic=cols[order[t]];
						lastIdx=t;
					}else if(!lastTraffic.equals(cols[order[t]])){
						for(int i=lastIdx;i<t;i++){
							rank[order[i]-2]=(double)(lastIdx+1+t)/2;
						}
						lastTraffic=cols[order[t]];
						lastIdx=t;
					}
					if(t==23){
						for(int i=lastIdx;i<=t;i++){
							rank[order[i]-2]=(double)(lastIdx+2+t)/2;
						}
					}
				}
				for(int t=0;t<24;t++){
					rank[t]=(rank[t]-12.5)/(Math.sqrt(Math.pow(23, 2)/12));
				}
					
				trafficRank.put(cols[0], rank);
				
//				for(int t=0;t<24;t++) System.out.print(rank[t]+"  ");
//				System.out.println();System.out.println();
			}
			System.out.println(instances.size());
			ins.close();
			resultWriter.close();
		}catch(Exception e){
			e.printStackTrace();
		}
	}
	/*
	 * district:2884
	 * city:38520
	% maxTraffic: 24780023840
	% minTraffic: 64
	@ATTRIBUTE id STRING
	@ATTRIBUTE hourEntropy NUMERIC, district(min0, max4.481, mean3.603, stdDev0.656),city(min0, max4.559, mean3.492, stdDev0.802)
	@ATTRIBUTE traffic1 NUMERIC
	@ATTRIBUTE traffic2 NUMERIC
	@ATTRIBUTE traffic3 NUMERIC
	@ATTRIBUTE traffic4 NUMERIC
	@ATTRIBUTE traffic5 NUMERIC
	@ATTRIBUTE traffic6 NUMERIC
	@ATTRIBUTE traffic7 NUMERIC
	@ATTRIBUTE traffic8 NUMERIC
	@ATTRIBUTE traffic9 NUMERIC
	@ATTRIBUTE traffic10 NUMERIC
	@ATTRIBUTE traffic11 NUMERIC
	@ATTRIBUTE traffic12 NUMERIC
	@ATTRIBUTE traffic13 NUMERIC
	@ATTRIBUTE traffic14 NUMERIC
	@ATTRIBUTE traffic15 NUMERIC
	@ATTRIBUTE traffic16 NUMERIC
	@ATTRIBUTE traffic17 NUMERIC
	@ATTRIBUTE traffic18 NUMERIC
	@ATTRIBUTE traffic19 NUMERIC
	@ATTRIBUTE traffic20 NUMERIC
	@ATTRIBUTE traffic21 NUMERIC
	@ATTRIBUTE traffic22 NUMERIC
	@ATTRIBUTE traffic23 NUMERIC
	@ATTRIBUTE traffic24 NUMERIC
	@ATTRIBUTE userEntropy NUMERIC, district(min0,max8.59,mean4.275,stdDev1.6), city(min0,max9.266,mean3.694,stdDev1.629)
	@ATTRIBUTE typeEntropy NUMERIC, district(min0,max2.166,mean1.061,stdDev0.366), city(min0,max2.391,mean1.08,stdDev0.378)
	@ATTRIBUTE trafficSum NUMERIC, district(min1368,max24780023840,mean525449103,stdDev949833404), city(min64,max24780023840,mean305072764,stdDev606101636)
	 */
	public static double distance(String[] first,String[] second){
		double d2=0D;
		d2+=Math.pow((Double.parseDouble(first[1])/4.56-Double.parseDouble(second[1])/4.56),2);   //hourEntropy
		d2+=Math.pow((Double.parseDouble(first[26])/9.27-Double.parseDouble(second[26])/9.27),2);  //userEntropy
		d2+=Math.pow((Double.parseDouble(first[27])/2.40-Double.parseDouble(second[27])/2.40),2);  //typeEntropy
		//d2+=Math.pow((Double.parseDouble(first[28])/2.48E10-Double.parseDouble(second[28])/2.48E10),2);  //sum traffic
		//d2+=(1-pearsonRho(trafficRank.get(first[0]),trafficRank.get(second[0])))/2;   //times correlation, Spearman's r
		d2+=(1-pearsonRho(normalizedTraffic.get(first[0]),normalizedTraffic.get(second[0])))/2;   //times correlation, pearson's r
		/*
		 * for pearsonRho (only for standardized data)
		 * www.econ.upf.edu/~michael/stanford/maeb6.pdf
		 * 
		 * d=sqrt(2(1-r)),  0<=d<=sqrt(2)
		 * normalized:
		 * dnorm=sqrt((1-r)/2), 0<=dnorm<=1
		 */
		return Math.sqrt(d2);
	}
	/*
	 * Gaussian kernel
	 * a(i,j)=exp(-||xi-xj||/2σ^2)
	 */
	public static double affinity(double d,double sigma){
		double a=Math.exp(-Math.pow(d,2)/(2*Math.pow(sigma,2)));
		if(a<affinityThreshold)
			return 0;
		else 
			return a;
	}
	
	

	/*
	 * ρX,Y=cov(X,Y)/(σX*σY)=E[(X-u)(Y-v)]/(σX*σY)
	 * ρ=Σ(Xi-u)(Yi-v)/(Math.sqrt(Σ(Xi-u)^2)*Math.sqrt(Σ(Yi-u)^2))
	 */
	public static Double pearsonRho(double[] first, double[] second){
		double u = 0,v = 0;
		for(int t=0;t<first.length;t++){
			u+=first[t];
			v+=second[t];
		}
		u=u/first.length;
		v=v/first.length;
		double e=0,vx=0,vy=0;
		for(int t=0;t<first.length;t++){
			e+=(first[t]-u)*(second[t]-v);
			vx+=Math.pow((first[t]-u),2);
			vy+=Math.pow((second[t]-v),2);
		}
		return e/(Math.sqrt(vx)*Math.sqrt(vy));
	}
	public static double spearmanRho(double[] firstRank,double[] secondRank){
		// or return pearsonRho(firstRank,secondRank);
		double sigmaD2=0D;
		for(int t=0;t<firstRank.length;t++)
			sigmaD2+=Math.pow(firstRank[t]-secondRank[t],2);
		return 1-(6*sigmaD2)/(firstRank.length*(Math.pow(firstRank.length, 2)-1));
	}
	
	
	

	//just for test
	@SuppressWarnings("unused")
	private static void loadInstancesTest(){
		for(int t=0;t<6;t++){
			instances.put(String.valueOf(t), new String[]{String.valueOf(t)});
			instanceKeyList.add(String.valueOf(t));
		}
	}
	@SuppressWarnings("unused")
	private static Double affinityTest(String[] first,String[] second){
		if(first[0].equals("0")&&second[0].equals("1"))return 0.8D;
		if(first[0].equals("0")&&second[0].equals("2"))return 0.6D;
		if(first[0].equals("0")&&second[0].equals("4"))return 0.1D;
		if(first[0].equals("1")&&second[0].equals("0"))return 0.8D;
		if(first[0].equals("1")&&second[0].equals("2"))return 0.8D;
		if(first[0].equals("2")&&second[0].equals("0"))return 0.6D;
		if(first[0].equals("2")&&second[0].equals("1"))return 0.8D;
		if(first[0].equals("2")&&second[0].equals("3"))return 0.2D;
		if(first[0].equals("3")&&second[0].equals("2"))return 0.2D;
		if(first[0].equals("3")&&second[0].equals("4"))return 0.8D;
		if(first[0].equals("3")&&second[0].equals("5"))return 0.7D;
		if(first[0].equals("4")&&second[0].equals("0"))return 0.1D;
		if(first[0].equals("4")&&second[0].equals("3"))return 0.8D;
		if(first[0].equals("4")&&second[0].equals("5"))return 0.8D;
		if(first[0].equals("5")&&second[0].equals("3"))return 0.7D;
		if(first[0].equals("5")&&second[0].equals("4"))return 0.8D;
		return 0D;
	}
	@SuppressWarnings("unused")
	private static void testDistance(){
		/*
		 * mean: 0.4840443715507323
		 * variance: 0.029066320442783875
		 * max: 1.5596958357717703
		 * min: 0.050596193951451034
		 */
		double sum=0,sum2=0,dis,count=0;
		double uu=0,max=0,min=2;

		FileWriter resultWriter;
		try {
			resultWriter = new FileWriter(getFilename()+".distance");
			
			for(int i=0;i<instanceKeyList.size();i++){
				for(int j=0;j<i;j++){
					dis=distance(instances.get(instanceKeyList.get(i)),instances.get(instanceKeyList.get(j)));
					sum+=dis;
					sum2+=Math.pow(dis, 2);
					if(max<dis)max=dis;
					if(min>dis)min=dis;
					count++;
					resultWriter.write(dis+"\r\n");
				}
			}
			double u=sum/count;
			uu=sum2/count-Math.pow(u, 2);
			System.out.println("mean: "+u);
			System.out.println("variance: "+uu);
			System.out.println("deviation: "+Math.sqrt(uu));
			System.out.println("max: "+max);
			System.out.println("min: "+min);
			System.out.println("count: "+count);
			resultWriter.close();
			resultWriter = new FileWriter(getFilename()+".distance.summary");
			resultWriter.write("%mean: "+u+"\r\n");
			resultWriter.write("%variance: "+uu+"\r\n");
			resultWriter.write("%deviation: "+Math.sqrt(uu)+"\r\n");
			resultWriter.write("%max: "+max+"\r\n");
			resultWriter.write("%min: "+min+"\r\n");
			resultWriter.write("%count: "+count+"\r\n");
			resultWriter.close();
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
	}
}

/**
 * Reference
 *
 * Spectral clustering method
[Shi and Malik 2000, Scott and Longuet-Higgins, Ng et al. 2002]
? Use k smallest eigenvectors
? Directly compute k-way partitioning
? Usually performs better


Relaxing RatioCut leads to unormalised spectral clustering [von Luxburg 2007]
Compute Laplacian L = D-A  (unormalised)

**/